/*
 * SteinhartHartModel.cpp
 *
 * Copyright 2012 Pieter Agten
 *
 * This file is part of tltgen.
 *
 * tltgen is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * tltgen is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with tltgen.  If not, see <http://www.gnu.org/licenses/>.
 */

/**
 * @file   SteinhartHartModel.cpp
 * @author Pieter Agten (pieter.agten@gmail.com)
 * @date   2 sep 2012
 */

#include <cmath>
#include <stdexcept>
#include <iostream>
#include <list>
#include <boost/lexical_cast.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/make_shared.hpp>

#include "SteinhartHartModelThermistor.hpp"
#include "LeastSquaresSolver.hpp"
#include "ThermistorTable.hpp"
#include "ThermistorTableRow.hpp"

#include "SteinhartHartModel.hpp"

#define DEBUG false

using boost::shared_ptr;
using boost::make_shared;

namespace tltgen {

  SteinhartHartModel::SteinhartHartModel() { }

  SteinhartHartModel::~SteinhartHartModel() { }

  std::string SteinhartHartModel::getName() const
  {
    return "Steinhart-Hart model";
  }
  unsigned int SteinhartHartModel::getMinNbValues() const
  {
    return 3;
  }

  shared_ptr<Thermistor> SteinhartHartModel::generate(const shared_ptr<ThermistorTable>& thermistorTable, const LeastSquaresSolver& solver)
  {
    double params[3];
    solver.calculateParameters(this, thermistorTable.get(), params);

    double a = params[0];
    double b = params[1];
    double c = params[2];

    return make_shared<SteinhartHartModelThermistor>(a, b, c);
  }

  unsigned int SteinhartHartModel::getNbParameters() const {
    return 3;
  }

  std::string SteinhartHartModel::getParameterName(unsigned int pi) const
  {
    if (! isValidParameterIndex(pi)) {
      throw std::invalid_argument("Invalid parameter index specified.");
    }

    switch (pi) {
    case 0:
      return "A";
    case 1:
      return "B";
    case 2:
      return "C";
    default:
      throw std::runtime_error("Default switch case should not be reachable.");
    }
  }

  ThermistorModel::ParameterType SteinhartHartModel::getParameterType(unsigned int pi) const
  {
    if (! isValidParameterIndex(pi)) {
      throw std::invalid_argument("Invalid parameter index specified.");
    }

    switch (pi) {
    case 0:
      return ThermistorModel::OTHER;
    case 1:
      return ThermistorModel::OTHER;
    case 2:
      return ThermistorModel::OTHER;
    default:
      throw std::runtime_error("Default switch case should not be reachable.");
    }
  }

  void SteinhartHartModel::getInitialGuess(double *guess) const
  {
    guess[0] = 1E-3;
    guess[1] = 1E-4;
    guess[2] = 1E-7;
  }
    
  bool SteinhartHartModel::isParameterFixed(unsigned int pi) const
  {
    if (! isValidParameterIndex(pi)) {
      throw std::invalid_argument(boost::lexical_cast<std::string>(pi) + " is not a valid parameter index for this Steinhart-Hart model generator.");
    }

    return false;
  }

  bool SteinhartHartModel::hasParameterLowerBound(unsigned int pi) const
  {
    if (! isValidParameterIndex(pi)) {
      throw std::invalid_argument(boost::lexical_cast<std::string>(pi) + " is not a valid parameter index for this Steinhart-Hart model generator.");
    }

    return true;
  }

  double SteinhartHartModel::getParameterLowerBound(unsigned int pi) const
  {
    if (! isValidParameterIndex(pi)) {
      throw std::invalid_argument(boost::lexical_cast<std::string>(pi) + " is not a valid parameter index for this Steinhart-Hart model generator.");
    }

    switch (pi) {
    case 0:
      // Minimum A value
      return 1E-14;
    case 1:
      // Minimum B value
      return 1E-14;
    case 2:
      // Minimum C value
      return 1E-14;
    default:
      throw std::runtime_error("Default switch case should not be reachable.");
    }
  }

  bool SteinhartHartModel::hasParameterUpperBound(unsigned int pi) const
  {
    if (! isValidParameterIndex(pi)) {
      throw std::invalid_argument(boost::lexical_cast<std::string>(pi) + " is not a valid parameter index for this Steinhart-Hart model generator.");
    }

    return false;
  }
  double SteinhartHartModel::getParameterUpperBound(unsigned int pi) const
  {
    throw std::runtime_error("This Steinhart-Hart model generator has no upper bound on parameter " + boost::lexical_cast<std::string>(pi) + ".");
  }

  void SteinhartHartModel::applyModelFunction(const ThermistorTable* thermistorTable, const double *params, double *residues) const
  {
    double a = params[0];
    double b = params[1];
    double c = params[2];

    SteinhartHartModelThermistor thermistor(a, b, c);

#if DEBUG
    std::cout << std::endl << "Trying parameters " << a << ", " << b << ", " << c << "." << std::endl;
#endif
    unsigned int i = 0;
    for (std::list<ThermistorTableRow *>::const_iterator it = thermistorTable->begin(); it != thermistorTable->end(); ++it) {
      if ((*it)->isValid()) {
	double Ri = (*it)->getResistance();
	double Ti = thermistor.getTemperature(Ri);
	double residue = Ti - (*it)->getTemperature();
#if DEBUG
	std::cout << "i=" << i << " R=" << Ri << " Ti=" << (*it)->getTemperature() << " Ti'=" << Ti << " Residue=" << residue << std::endl;
#endif
	residues[i] = residue;
	i += 1;
      }
    }
  }


  /**  
   * Jacobian matrix J(i,j) = dfi/dxj
   *   where fi = (Ti - ti)
   *         Ti = 1.0/(A + B*ln(ri) + C*ln³(ri))
   *   the xj are the parameters (A,B,C)
   *
   *  dfi/dA = -Ti²     dfi/dB = -ln(ri)*Ti²     dfi/dB = -ln³(ri)*Ti²
   */
  void SteinhartHartModel::applyModelDerivative(const ThermistorTable* thermistorTable, const double *params, double *jacobian) const
  {
    double a = params[0];
    double b = params[1];
    double c = params[2];

    SteinhartHartModelThermistor thermistor(a, b, c);

    unsigned int i = 0;
    for (std::list<ThermistorTableRow *>::const_iterator it = thermistorTable->begin(); it != thermistorTable->end(); ++it) {
      if ((*it)->isValid()) {
	double Ri = (*it)->getResistance();
	double ln_ri = log(Ri);
	double Ti2 = pow(thermistor.getTemperature(Ri), 2);
    
	jacobian[i * 3 + 0] = -Ti2;
	jacobian[i * 3 + 1] = -ln_ri * Ti2;
	jacobian[i * 3 + 2] = -pow(ln_ri, 3) * Ti2;
	i += 1;
      }
    }
  }

}

